Sensory Evaluation of French-fries

### Project Description The client wants to examine the effects of two factors (types of oil; age of oil) on the taste of French-fries. Her theory is that fries which were cooked in oil that has been reused over several days taste better than those cooked in fresh oil. She intends to test the effect of different types and ages of oil on the taste of the fries. She is also interested in whether certain features, such as texture, temperature, and color, are related to the perceived quality of the fries.

Her experiment has already been performed with four types of oil and each ranging in age from one to five days. The experiment ran for two weeks. At lunchtime each day, students came into the lab and were randomly given two boxes of fries. During week 1, the students were given one box with oil A and one box with oil B. During week 2, they were given one box with oil C and one box with oil D. After tasting the fries, students filled out a questionnaire in which they scored various attributes on a scale of 1 to 9. The attributes were: appearance, temperature, color, taste, texture, and overall liking. The questionnaire also recorded which fries they preferred, along with health information, allergies, and gender. They gathered a large amount of data, roughly 300-400 responses per day over the two weeks. The data were entered in Microsoft Excel and client agreed to analyze the data using whatever statistical software the consultants prefer.

Data Preparation

library(readxl)
library(plyr)
library(nlme)
library(dplyr)

## 
## Attaching package: 'dplyr'
## 
## The following object is masked from 'package:nlme':
## 
##     collapse
## 
## The following objects are masked from 'package:plyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## 
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## 
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

# function that reads all sheets in the Excel file and puts them in a list
read_excel_allsheets <- function(filename) {
  sheets <- readxl::excel_sheets(filename)
  x <- lapply(sheets, function(x) readxl::read_excel(filename, sheet = x))
  names(x) <- sheets
  x
}

# collapse the list of sheets into a single data frame
data <- ldply(read_excel_allsheets('data.xlsx'))
# fix column names
names(data) <- gsub('\\-|\\#| |\\/','.',names(data))
# select columns we're interested in
data <- select(data, Samp.Desc,Age,Temperature,Appearance,Color,Taste,Texture,Overall.Liking,Samp.Set)
data$Age <- as.factor(data$Age) # change Age from numeric to factor
names(data)[1]<-"Type"; names(data)[8]<-"Liking" # rename Samp.Desc to Type
names(data)[9] <- "Subject"
data$Subject <- rep(1:(length(data$Subject)/2), times=rep(2, length(data$Subject)/2))


str(data) # check the data

## 'data.frame':    1890 obs. of  9 variables:
##  $ Type       : chr  "Mel Fry Free" "Clear Valley" "Clear Valley" "Mel Fry Free" ...
##  $ Age        : Factor w/ 5 levels "1","2","3","4",..: 2 2 2 2 4 4 2 2 5 5 ...
##  $ Temperature: num  5 5 5 6 5 4 5 4 5 4 ...
##  $ Appearance : num  3 3 7 7 8 4 4 3 8 5 ...
##  $ Color      : num  3 2 5 5 5 6 4 5 5 5 ...
##  $ Taste      : num  3 4 7 7 8 6 5 4 2 2 ...
##  $ Texture    : num  3 3 7 7 7 6 7 7 5 5 ...
##  $ Liking     : num  5 4 7 7 8 6 6 4 2 2 ...
##  $ Subject    : int  1 1 2 2 3 3 4 4 5 5 ...

# assign new values to Temperature so that the higher the number, the more likeness; 
# 9 changed to 1, 8 changed to 2, 7 changed to 3, 6 changed to 4; ONE problem is the modified data are skewed extremely; I am struggling to find a way to solve this; any idea? the same issue happens to Color also
data$Temperature[data$Temperature==9] <- 1; data$Temperature[data$Temperature==8] <- 2;
data$Temperature[data$Temperature==7] <- 3; data$Temperature[data$Temperature==6] <- 4;
# assign new values to Color so that the higher the number, the more likeness; 
# 9 changed to 1, 8 changed to 2, 7 changed to 3, 6 changed to 4
data$Color[data$Color==9] <- 1; data$Color[data$Color==8] <- 2;
data$Color[data$Color==7] <- 3; data$Color[data$Color==6] <- 4;

Exploratory Data Analysis

# The histograms below show that the four variables (Appearance, Taste, Texture and Temperature) are approximately normal distributed.
par(mfrow=c(2,2), mar=c(4, 4, 2, 1))
hist(data$Appearance, col="green"); hist(data$Taste, col="green")
hist(data$Texture, col="green"); hist(data$Liking, col="green")

# data skewed extremely; I am struggling to find a way to solve this; any idea?
hist(data$Temperature, col="green"); hist(data$Color, col="green") 
# Scatterplots of every combination of attributes (to check for collinearity
pairs(~Temperature+Appearance+Color+Taste+Texture+Liking,data, main="Simple Scatterplot Matrix")

Question 1: How do oil types and ages affect the overall liking of French fries?

############ Overall Liking ############
# repeated measures: take into account the variation in taste between students (since each student tasted two boxes of fries)
repmeas <- lme(Liking ~ Type * Age, data=data, random = ~ 1 | Subject, method='ML'); summary(repmeas)

## Linear mixed-effects model fit by maximum likelihood
##  Data: data 
##        AIC      BIC    logLik
##   6400.395 6522.371 -3178.198
## 
## Random effects:
##  Formula: ~1 | Subject
##         (Intercept)  Residual
## StdDev:    0.985983 0.9891079
## 
## Fixed effects: Liking ~ Type * Age 
##                           Value Std.Error  DF  t-value p-value
## (Intercept)            6.093023 0.0957554 940 63.63109  0.0000
## TypeClear Valley      -0.116461 0.1567493 930 -0.74297  0.4577
## TypeMel Fry           -0.037209 0.0959068 930 -0.38797  0.6981
## TypeMel Fry Free      -0.038336 0.1567493 930 -0.24457  0.8068
## Age2                   0.086977 0.1493703 940  0.58229  0.5605
## Age3                   0.078405 0.2559172 940  0.30637  0.7594
## Age4                   0.281977 0.1911781 940  1.47494  0.1406
## Age5                   0.906977 0.2153243 940  4.21214  0.0000
## TypeClear Valley:Age2  0.088842 0.2376711 930  0.37380  0.7086
## TypeMel Fry:Age2       0.163876 0.1496065 930  1.09538  0.2736
## TypeMel Fry Free:Age2  0.182145 0.2376711 930  0.76638  0.4436
## TypeClear Valley:Age3 -0.138301 0.3683134 930 -0.37550  0.7074
## TypeMel Fry:Age3       0.180066 0.2563218 930  0.70250  0.4825
## TypeMel Fry Free:Age3 -0.021982 0.3683134 930 -0.05968  0.9524
## TypeClear Valley:Age4  0.465145 0.2790864 930  1.66667  0.0959
## TypeMel Fry:Age4       0.078876 0.1914804 930  0.41193  0.6805
## TypeMel Fry Free:Age4  0.123862 0.2790864 930  0.44381  0.6573
## TypeClear Valley:Age5 -0.003539 0.2967331 930 -0.01193  0.9905
## TypeMel Fry:Age5       0.112681 0.2156647 930  0.52248  0.6015
## TypeMel Fry Free:Age5 -0.268331 0.2967331 930 -0.90428  0.3661
##  Correlation: 
##                       (Intr) TypClV TypMlF TypMFF Age2   Age3   Age4  
## TypeClear Valley      -0.611                                          
## TypeMel Fry           -0.501  0.306                                   
## TypeMel Fry Free      -0.611  0.686  0.306                            
## Age2                  -0.641  0.392  0.321  0.392                     
## Age3                  -0.374  0.229  0.187  0.229  0.240              
## Age4                  -0.501  0.306  0.251  0.306  0.321  0.187       
## Age5                  -0.445  0.272  0.223  0.272  0.285  0.166  0.223
## TypeClear Valley:Age2  0.403 -0.660 -0.202 -0.452 -0.628 -0.151 -0.202
## TypeMel Fry:Age2       0.321 -0.196 -0.641 -0.196 -0.501 -0.120 -0.161
## TypeMel Fry Free:Age2  0.403 -0.452 -0.202 -0.660 -0.628 -0.151 -0.202
## TypeClear Valley:Age3  0.260 -0.426 -0.130 -0.292 -0.167 -0.695 -0.130
## TypeMel Fry:Age3       0.187 -0.114 -0.374 -0.114 -0.120 -0.501 -0.094
## TypeMel Fry Free:Age3  0.260 -0.292 -0.130 -0.426 -0.167 -0.695 -0.130
## TypeClear Valley:Age4  0.343 -0.562 -0.172 -0.385 -0.220 -0.128 -0.685
## TypeMel Fry:Age4       0.251 -0.153 -0.501 -0.153 -0.161 -0.094 -0.501
## TypeMel Fry Free:Age4  0.343 -0.385 -0.172 -0.562 -0.220 -0.128 -0.685
## TypeClear Valley:Age5  0.323 -0.528 -0.162 -0.362 -0.207 -0.121 -0.162
## TypeMel Fry:Age5       0.223 -0.136 -0.445 -0.136 -0.143 -0.083 -0.112
## TypeMel Fry Free:Age5  0.323 -0.362 -0.162 -0.528 -0.207 -0.121 -0.162
##                       Age5   TCV:A2 TMF:A2 TMFF:A2 TCV:A3 TMF:A3 TMFF:A3
## TypeClear Valley                                                        
## TypeMel Fry                                                             
## TypeMel Fry Free                                                        
## Age2                                                                    
## Age3                                                                    
## Age4                                                                    
## Age5                                                                    
## TypeClear Valley:Age2 -0.179                                            
## TypeMel Fry:Age2      -0.143  0.315                                     
## TypeMel Fry Free:Age2 -0.179  0.697  0.315                              
## TypeClear Valley:Age3 -0.116  0.281  0.083  0.192                       
## TypeMel Fry:Age3      -0.083  0.075  0.240  0.075   0.348               
## TypeMel Fry Free:Age3 -0.116  0.192  0.083  0.281   0.741  0.348        
## TypeClear Valley:Age4 -0.153  0.370  0.110  0.254   0.239  0.064  0.164 
## TypeMel Fry:Age4      -0.112  0.101  0.321  0.101   0.065  0.187  0.065 
## TypeMel Fry Free:Age4 -0.153  0.254  0.110  0.370   0.164  0.064  0.239 
## TypeClear Valley:Age5 -0.726  0.348  0.104  0.239   0.225  0.060  0.154 
## TypeMel Fry:Age5      -0.501  0.090  0.285  0.090   0.058  0.166  0.058 
## TypeMel Fry Free:Age5 -0.726  0.239  0.104  0.348   0.154  0.060  0.225 
##                       TCV:A4 TMF:A4 TMFF:A4 TCV:A5 TMF:A5
## TypeClear Valley                                         
## TypeMel Fry                                              
## TypeMel Fry Free                                         
## Age2                                                     
## Age3                                                     
## Age4                                                     
## Age5                                                     
## TypeClear Valley:Age2                                    
## TypeMel Fry:Age2                                         
## TypeMel Fry Free:Age2                                    
## TypeClear Valley:Age3                                    
## TypeMel Fry:Age3                                         
## TypeMel Fry Free:Age3                                    
## TypeClear Valley:Age4                                    
## TypeMel Fry:Age4       0.343                             
## TypeMel Fry Free:Age4  0.734  0.343                      
## TypeClear Valley:Age5  0.297  0.081  0.203               
## TypeMel Fry:Age5       0.076  0.223  0.076   0.363       
## TypeMel Fry Free:Age5  0.203  0.081  0.297   0.763  0.363
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.04154704 -0.44142058  0.00821523  0.56182165  2.69036815 
## 
## Number of Observations: 1890
## Number of Groups: 945

# nothing significant in the interaction terms, so we probably can leave them out
repmeas2 <- lme(Liking ~ Type + Age, data=data, random = ~ 1 | Subject, method='ML'); summary(repmeas2)

## Linear mixed-effects model fit by maximum likelihood
##  Data: data 
##        AIC      BIC    logLik
##   6387.365 6442.808 -3183.683
## 
## Random effects:
##  Formula: ~1 | Subject
##         (Intercept)  Residual
## StdDev:   0.9860901 0.9933616
## 
## Fixed effects: Liking ~ Type + Age 
##                      Value  Std.Error  DF  t-value p-value
## (Intercept)       6.046374 0.07832796 942 77.19306  0.0000
## TypeClear Valley -0.021627 0.09266087 942 -0.23340  0.8155
## TypeMel Fry       0.043810 0.06144178 942  0.71302  0.4760
## TypeMel Fry Free -0.019246 0.09266087 942 -0.20771  0.8355
## Age2              0.190332 0.10036750 940  1.89635  0.0582
## Age3              0.079950 0.15855454 940  0.50425  0.6142
## Age4              0.450086 0.11985261 940  3.75533  0.0002
## Age5              0.847155 0.12682253 940  6.67985  0.0000
##  Correlation: 
##                  (Intr) TypClV TypMlF TypMFF Age2   Age3   Age4  
## TypeClear Valley -0.462                                          
## TypeMel Fry      -0.392  0.332                                   
## TypeMel Fry Free -0.462  0.725  0.332                            
## Age2             -0.534 -0.027  0.000 -0.027                     
## Age3             -0.320 -0.059  0.000 -0.059  0.272              
## Age4             -0.421 -0.082  0.000 -0.082  0.360  0.232       
## Age5             -0.380 -0.117  0.000 -0.117  0.341  0.223  0.295
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -3.097076824 -0.410551402  0.003857149  0.593343904  2.640946939 
## 
## Number of Observations: 1890
## Number of Groups: 945

# looks like there is no significant difference between oil types. Also no significant difference between ages 1, 2, and 3. Ages 4 and 5 are different, however. This suggests that type doesn't have an effect on liking, while age does. Specifically, oil that is > 3 days old corresponds with an increase in liking.
intervals(repmeas2)

## Approximate 95% confidence intervals
## 
##  Fixed effects:
##                         lower        est.     upper
## (Intercept)       5.892982504  6.04637431 6.1997661
## TypeClear Valley -0.203087648 -0.02162729 0.1598331
## TypeMel Fry      -0.076513624  0.04380952 0.1641327
## TypeMel Fry Free -0.200706695 -0.01924634 0.1622140
## Age2             -0.006220822  0.19033217 0.3868852
## Age3             -0.230552137  0.07995046 0.3904531
## Age4              0.215374186  0.45008552 0.6847969
## Age5              0.598794691  0.84715545 1.0955162
## attr(,"label")
## [1] "Fixed effects:"
## 
##  Random Effects:
##   Level: Subject 
##                     lower      est.   upper
## sd((Intercept)) 0.9049461 0.9860901 1.07451
## 
##  Within-group standard error:
##     lower      est.     upper 
## 0.9441392 0.9933616 1.0451503

# 95% confidence intervals reflect this as well.
anova(repmeas, repmeas2)

##          Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## repmeas      1 22 6400.395 6522.371 -3178.198                        
## repmeas2     2 10 6387.365 6442.808 -3183.683 1 vs 2 10.96971  0.5315

# by a likelihood ratio test (p=0.53), leaving out the interaction terms didn't affect fit

# make some plots
boxplot(Liking ~ Type, data=data) # again we don't see much difference among different types of oil

boxplot(Liking ~ Age, data=data) # apparently we can see Age 5 is significantly higher than other ages

# make a plot showing how overall likings change with ages of oil grouped by different types. We can see generally there is an uptrend with the increase of oil ages within each type of oil.
data.1 <- aggregate(data[,3:8], by=list(data$Type, data$Age), mean) 
library(ggplot2)
d <- ggplot(data=data.1, aes(x=Group.2, y=Liking)) 
d <- d + geom_bar(fill="grey", color="grey",stat = "identity") +
  guides(fill=F) + 
  ggtitle('Mean of Overall Likings in Different Types and Ages of Oils') + 
  ylab('Liking') + 
  xlab('Age') + 
  theme(legend.position='none') + 
  facet_grid(. ~ Group.1) + 
  geom_text(aes(label=round(Liking, 0), size=1, hjust=0.5, vjust=-1))
print(d)

Question 2: How do oil types and ages affect the characteristics (temperature, appearance, color, taste and texture) of French fries?

We could repeat the steps in question 1 for the characteristics. (I leave Temperature and Color undone here; may need help)

########### Appearance ############
fit.Appearance <- lme(Appearance ~ Type*Age, data=data, random = ~ 1 | Subject, method='ML'); summary(fit.Appearance)

## Linear mixed-effects model fit by maximum likelihood
##  Data: data 
##        AIC      BIC    logLik
##   6242.664 6364.639 -3099.332
## 
## Random effects:
##  Formula: ~1 | Subject
##         (Intercept) Residual
## StdDev:   0.9848551  0.92911
## 
## Fixed effects: Appearance ~ Type * Age 
##                           Value Std.Error  DF  t-value p-value
## (Intercept)            6.260465 0.0928312 940 67.43922  0.0000
## TypeClear Valley       0.052035 0.1519624 930  0.34242  0.7321
## TypeMel Fry           -0.162791 0.0900893 930 -1.80699  0.0711
## TypeMel Fry Free       0.130160 0.1519624 930  0.85653  0.3919
## Age2                  -0.007132 0.1448088 940 -0.04925  0.9607
## Age3                   0.396678 0.2481019 940  1.59885  0.1102
## Age4                   0.392313 0.1853398 940  2.11672  0.0345
## Age5                   0.984818 0.2087486 940  4.71772  0.0000
## TypeClear Valley:Age2 -0.086321 0.2304130 930 -0.37463  0.7080
## TypeMel Fry:Age2       0.302791 0.1405315 930  2.15461  0.0314
## TypeMel Fry Free:Age2 -0.050160 0.2304130 930 -0.21770  0.8277
## TypeClear Valley:Age3 -0.403622 0.3570657 930 -1.13039  0.2586
## TypeMel Fry:Age3      -0.122924 0.2407737 930 -0.51054  0.6098
## TypeMel Fry Free:Age3 -0.509525 0.3570657 930 -1.42698  0.1539
## TypeClear Valley:Age4 -0.033760 0.2705635 930 -0.12478  0.9007
## TypeMel Fry:Age4       0.135013 0.1798654 930  0.75063  0.4531
## TypeMel Fry Free:Age4 -0.203990 0.2705635 930 -0.75395  0.4511
## TypeClear Valley:Age5 -0.057318 0.2876713 930 -0.19925  0.8421
## TypeMel Fry:Age5       0.049583 0.2025828 930  0.24476  0.8067
## TypeMel Fry Free:Age5 -0.348776 0.2876713 930 -1.21241  0.2257
##  Correlation: 
##                       (Intr) TypClV TypMlF TypMFF Age2   Age3   Age4  
## TypeClear Valley      -0.611                                          
## TypeMel Fry           -0.485  0.296                                   
## TypeMel Fry Free      -0.611  0.705  0.296                            
## Age2                  -0.641  0.392  0.311  0.392                     
## Age3                  -0.374  0.229  0.182  0.229  0.240              
## Age4                  -0.501  0.306  0.243  0.306  0.321  0.187       
## Age5                  -0.445  0.272  0.216  0.272  0.285  0.166  0.223
## TypeClear Valley:Age2  0.403 -0.660 -0.195 -0.465 -0.628 -0.151 -0.202
## TypeMel Fry:Age2       0.311 -0.190 -0.641 -0.190 -0.485 -0.116 -0.156
## TypeMel Fry Free:Age2  0.403 -0.465 -0.195 -0.660 -0.628 -0.151 -0.202
## TypeClear Valley:Age3  0.260 -0.426 -0.126 -0.300 -0.167 -0.695 -0.130
## TypeMel Fry:Age3       0.182 -0.111 -0.374 -0.111 -0.116 -0.485 -0.091
## TypeMel Fry Free:Age3  0.260 -0.300 -0.126 -0.426 -0.167 -0.695 -0.130
## TypeClear Valley:Age4  0.343 -0.562 -0.166 -0.396 -0.220 -0.128 -0.685
## TypeMel Fry:Age4       0.243 -0.148 -0.501 -0.148 -0.156 -0.091 -0.485
## TypeMel Fry Free:Age4  0.343 -0.396 -0.166 -0.562 -0.220 -0.128 -0.685
## TypeClear Valley:Age5  0.323 -0.528 -0.157 -0.372 -0.207 -0.121 -0.162
## TypeMel Fry:Age5       0.216 -0.132 -0.445 -0.132 -0.138 -0.081 -0.108
## TypeMel Fry Free:Age5  0.323 -0.372 -0.157 -0.528 -0.207 -0.121 -0.162
##                       Age5   TCV:A2 TMF:A2 TMFF:A2 TCV:A3 TMF:A3 TMFF:A3
## TypeClear Valley                                                        
## TypeMel Fry                                                             
## TypeMel Fry Free                                                        
## Age2                                                                    
## Age3                                                                    
## Age4                                                                    
## Age5                                                                    
## TypeClear Valley:Age2 -0.179                                            
## TypeMel Fry:Age2      -0.138  0.305                                     
## TypeMel Fry Free:Age2 -0.179  0.715  0.305                              
## TypeClear Valley:Age3 -0.116  0.281  0.081  0.198                       
## TypeMel Fry:Age3      -0.081  0.073  0.240  0.073   0.337               
## TypeMel Fry Free:Age3 -0.116  0.198  0.081  0.281   0.756  0.337        
## TypeClear Valley:Age4 -0.153  0.370  0.107  0.261   0.239  0.062  0.168 
## TypeMel Fry:Age4      -0.108  0.098  0.321  0.098   0.063  0.187  0.063 
## TypeMel Fry Free:Age4 -0.153  0.261  0.107  0.370   0.168  0.062  0.239 
## TypeClear Valley:Age5 -0.726  0.348  0.100  0.246   0.225  0.059  0.158 
## TypeMel Fry:Age5      -0.485  0.087  0.285  0.087   0.056  0.166  0.056 
## TypeMel Fry Free:Age5 -0.726  0.246  0.100  0.348   0.158  0.059  0.225 
##                       TCV:A4 TMF:A4 TMFF:A4 TCV:A5 TMF:A5
## TypeClear Valley                                         
## TypeMel Fry                                              
## TypeMel Fry Free                                         
## Age2                                                     
## Age3                                                     
## Age4                                                     
## Age5                                                     
## TypeClear Valley:Age2                                    
## TypeMel Fry:Age2                                         
## TypeMel Fry Free:Age2                                    
## TypeClear Valley:Age3                                    
## TypeMel Fry:Age3                                         
## TypeMel Fry Free:Age3                                    
## TypeClear Valley:Age4                                    
## TypeMel Fry:Age4       0.332                             
## TypeMel Fry Free:Age4  0.750  0.332                      
## TypeClear Valley:Age5  0.297  0.078  0.209               
## TypeMel Fry:Age5       0.074  0.223  0.074   0.352       
## TypeMel Fry Free:Age5  0.209  0.078  0.297   0.777  0.352
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.53209930 -0.48450265  0.07061119  0.52055717  2.78199021 
## 
## Number of Observations: 1890
## Number of Groups: 945

# No significant interaction terms suggests we do not have to worry about interaction. So we can repeat the process but with a simple additive model.
fit.Appearance2 <- lme(Appearance ~ Type + Age, data=data, random = ~ 1 | Subject, method='ML'); summary(fit.Appearance2)

## Linear mixed-effects model fit by maximum likelihood
##  Data: data 
##        AIC      BIC    logLik
##   6230.606 6286.049 -3105.303
## 
## Random effects:
##  Formula: ~1 | Subject
##         (Intercept)  Residual
## StdDev:   0.9844683 0.9337987
## 
## Fixed effects: Appearance ~ Type + Age 
##                      Value  Std.Error  DF  t-value p-value
## (Intercept)       6.266142 0.07622915 942 82.20139  0.0000
## TypeClear Valley -0.007263 0.08984364 942 -0.08084  0.9356
## TypeMel Fry      -0.060952 0.05775768 942 -1.05531  0.2916
## TypeMel Fry Free -0.012025 0.08984364 942 -0.13385  0.8936
## Age2              0.059677 0.09827118 940  0.60727  0.5438
## Age3              0.155180 0.15524290 940  0.99959  0.3178
## Age4              0.385393 0.11734932 940  3.28415  0.0011
## Age5              0.908378 0.12417366 940  7.31538  0.0000
##  Correlation: 
##                  (Intr) TypClV TypMlF TypMFF Age2   Age3   Age4  
## TypeClear Valley -0.459                                          
## TypeMel Fry      -0.379  0.321                                   
## TypeMel Fry Free -0.459  0.742  0.321                            
## Age2             -0.537 -0.027  0.000 -0.027                     
## Age3             -0.321 -0.059  0.000 -0.059  0.272              
## Age4             -0.424 -0.082  0.000 -0.082  0.360  0.232       
## Age5             -0.383 -0.118  0.000 -0.118  0.341  0.223  0.295
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.59062546 -0.46304456  0.08888843  0.55693966  2.76946836 
## 
## Number of Observations: 1890
## Number of Groups: 945

# The above result shows the factor "Type" does not have a significant impact on the Appearance of French fries, but the factor "Age" does. Again, ages of 4 and 5 are significant.
# Use diagnostic plots to check heteroscedasticity, normality, and influential observerations. These plots appears to indicate the model is a good fit.
plot(fit.Appearance2); qqnorm(resid(fit.Appearance2)); qqline(resid(fit.Appearance2)) # diagnostic plots

anova(fit.Appearance, fit.Appearance2) # again (p=0.45) leaving out interactions was ok

##                 Model df      AIC      BIC    logLik   Test  L.Ratio
## fit.Appearance      1 22 6242.664 6364.639 -3099.332                
## fit.Appearance2     2 10 6230.606 6286.049 -3105.303 1 vs 2 11.94167
##                 p-value
## fit.Appearance         
## fit.Appearance2  0.4504

# make some plots
boxplot(Appearance ~ Type, data) # again we don't see much difference among different types of oil

boxplot(Appearance ~ Age, data) # apparently we can see Age 5 is significantly higher than other ages

# make a plot showing how Appearance change with ages of oil grouped by different types. We can see generally there is an uptrend with the increase of oil ages within each type of oil.
e <- ggplot(data=data.1, aes(x=Group.2, y=Appearance)) 
e <- e + geom_bar(fill="grey", color="grey",stat = "identity") +
  guides(fill=F) + 
  ggtitle('Mean of Appearance in Different Types and Ages of Oils') + 
  ylab('Appearance') + 
  xlab('Age') + 
  theme(legend.position='none') + 
  facet_grid(. ~ Group.1) + 
  geom_text(aes(label=round(Liking, 0), size=1, hjust=0.5, vjust=-1))
print(e)

########### Taste ############
fit.Taste <- lme(Taste ~ Type*Age, data=data, random = ~ 1 | Subject, method='ML'); summary(fit.Taste)

## Linear mixed-effects model fit by maximum likelihood
##  Data: data 
##        AIC      BIC    logLik
##   6492.408 6614.383 -3224.204
## 
## Random effects:
##  Formula: ~1 | Subject
##         (Intercept) Residual
## StdDev:   0.9991054 1.019066
## 
## Fixed effects: Taste ~ Type * Age 
##                           Value Std.Error  DF  t-value p-value
## (Intercept)            5.827907 0.0978487 940 59.56037  0.0000
## TypeClear Valley      -0.038844 0.1601759 930 -0.24251  0.8084
## TypeMel Fry            0.037209 0.0988117 930  0.37657  0.7066
## TypeMel Fry Free       0.093968 0.1601759 930  0.58666  0.5576
## Age2                   0.332093 0.1526357 940  2.17572  0.0298
## Age3                   0.429236 0.2615117 940  1.64136  0.1011
## Age4                   0.477649 0.1953574 940  2.44500  0.0147
## Age5                   1.247565 0.2200314 940  5.66994  0.0000
## TypeClear Valley:Age2 -0.025917 0.2428668 930 -0.10671  0.9150
## TypeMel Fry:Age2      -0.057209 0.1541378 930 -0.37116  0.7106
## TypeMel Fry Free:Age2 -0.044444 0.2428668 930 -0.18300  0.8548
## TypeClear Valley:Age3 -0.190521 0.3763650 930 -0.50621  0.6128
## TypeMel Fry:Age3       0.134219 0.2640853 930  0.50824  0.6114
## TypeMel Fry Free:Age3 -0.212222 0.3763650 930 -0.56387  0.5730
## TypeClear Valley:Age4  0.391184 0.2851874 930  1.37167  0.1705
## TypeMel Fry:Age4       0.073902 0.1972799 930  0.37460  0.7080
## TypeMel Fry Free:Age4 -0.083734 0.2851874 930 -0.29361  0.7691
## TypeClear Valley:Age5 -0.156627 0.3032199 930 -0.51655  0.6056
## TypeMel Fry:Age5      -0.037209 0.2221968 930 -0.16746  0.8670
## TypeMel Fry Free:Age5 -0.489440 0.3032199 930 -1.61414  0.1068
##  Correlation: 
##                       (Intr) TypClV TypMlF TypMFF Age2   Age3   Age4  
## TypeClear Valley      -0.611                                          
## TypeMel Fry           -0.505  0.308                                   
## TypeMel Fry Free      -0.611  0.680  0.308                            
## Age2                  -0.641  0.392  0.324  0.392                     
## Age3                  -0.374  0.229  0.189  0.229  0.240              
## Age4                  -0.501  0.306  0.253  0.306  0.321  0.187       
## Age5                  -0.445  0.272  0.225  0.272  0.285  0.166  0.223
## TypeClear Valley:Age2  0.403 -0.660 -0.203 -0.449 -0.628 -0.151 -0.202
## TypeMel Fry:Age2       0.324 -0.198 -0.641 -0.198 -0.505 -0.121 -0.162
## TypeMel Fry Free:Age2  0.403 -0.449 -0.203 -0.660 -0.628 -0.151 -0.202
## TypeClear Valley:Age3  0.260 -0.426 -0.131 -0.290 -0.167 -0.695 -0.130
## TypeMel Fry:Age3       0.189 -0.115 -0.374 -0.115 -0.121 -0.505 -0.095
## TypeMel Fry Free:Age3  0.260 -0.290 -0.131 -0.426 -0.167 -0.695 -0.130
## TypeClear Valley:Age4  0.343 -0.562 -0.173 -0.382 -0.220 -0.128 -0.685
## TypeMel Fry:Age4       0.253 -0.154 -0.501 -0.154 -0.162 -0.095 -0.505
## TypeMel Fry Free:Age4  0.343 -0.382 -0.173 -0.562 -0.220 -0.128 -0.685
## TypeClear Valley:Age5  0.323 -0.528 -0.163 -0.359 -0.207 -0.121 -0.162
## TypeMel Fry:Age5       0.225 -0.137 -0.445 -0.137 -0.144 -0.084 -0.112
## TypeMel Fry Free:Age5  0.323 -0.359 -0.163 -0.528 -0.207 -0.121 -0.162
##                       Age5   TCV:A2 TMF:A2 TMFF:A2 TCV:A3 TMF:A3 TMFF:A3
## TypeClear Valley                                                        
## TypeMel Fry                                                             
## TypeMel Fry Free                                                        
## Age2                                                                    
## Age3                                                                    
## Age4                                                                    
## Age5                                                                    
## TypeClear Valley:Age2 -0.179                                            
## TypeMel Fry:Age2      -0.144  0.317                                     
## TypeMel Fry Free:Age2 -0.179  0.692  0.317                              
## TypeClear Valley:Age3 -0.116  0.281  0.084  0.191                       
## TypeMel Fry:Age3      -0.084  0.076  0.240  0.076   0.351               
## TypeMel Fry Free:Age3 -0.116  0.191  0.084  0.281   0.736  0.351        
## TypeClear Valley:Age4 -0.153  0.370  0.111  0.252   0.239  0.065  0.163 
## TypeMel Fry:Age4      -0.112  0.102  0.321  0.102   0.066  0.187  0.066 
## TypeMel Fry Free:Age4 -0.153  0.252  0.111  0.370   0.163  0.065  0.239 
## TypeClear Valley:Age5 -0.726  0.348  0.104  0.237   0.225  0.061  0.153 
## TypeMel Fry:Age5      -0.505  0.090  0.285  0.090   0.058  0.166  0.058 
## TypeMel Fry Free:Age5 -0.726  0.237  0.104  0.348   0.153  0.061  0.225 
##                       TCV:A4 TMF:A4 TMFF:A4 TCV:A5 TMF:A5
## TypeClear Valley                                         
## TypeMel Fry                                              
## TypeMel Fry Free                                         
## Age2                                                     
## Age3                                                     
## Age4                                                     
## Age5                                                     
## TypeClear Valley:Age2                                    
## TypeMel Fry:Age2                                         
## TypeMel Fry Free:Age2                                    
## TypeClear Valley:Age3                                    
## TypeMel Fry:Age3                                         
## TypeMel Fry Free:Age3                                    
## TypeClear Valley:Age4                                    
## TypeMel Fry:Age4       0.346                             
## TypeMel Fry Free:Age4  0.729  0.346                      
## TypeClear Valley:Age5  0.297  0.082  0.202               
## TypeMel Fry:Age5       0.077  0.223  0.077   0.366       
## TypeMel Fry Free:Age5  0.202  0.082  0.297   0.759  0.366
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.58215055 -0.46914357  0.03328181  0.58533024  2.65439837 
## 
## Number of Observations: 1890
## Number of Groups: 945

# No significant interaction terms suggests we do not have to worry about interaction. So we can repeat the process but with a simple additive model.
fit.Taste2 <- lme(Taste ~ Type + Age, data=data, random = ~ 1 | Subject, method='ML'); summary(fit.Taste2)

## Linear mixed-effects model fit by maximum likelihood
##  Data: data 
##        AIC      BIC    logLik
##   6479.513 6534.956 -3229.756
## 
## Random effects:
##  Formula: ~1 | Subject
##         (Intercept) Residual
## StdDev:   0.9990248 1.023574
## 
## Fixed effects: Taste ~ Type + Age 
##                      Value  Std.Error  DF  t-value p-value
## (Intercept)       5.845304 0.07995253 942 73.10968  0.0000
## TypeClear Valley -0.008676 0.09468217 942 -0.09164  0.9270
## TypeMel Fry       0.036190 0.06331049 942  0.57163  0.5677
## TypeMel Fry Free -0.027724 0.09468217 942 -0.29281  0.7697
## Age2              0.302527 0.10227169 940  2.95807  0.0032
## Age3              0.366272 0.16156267 940  2.26706  0.0236
## Age4              0.580915 0.12212648 940  4.75667  0.0000
## Age5              1.060212 0.12922864 940  8.20415  0.0000
##  Correlation: 
##                  (Intr) TypClV TypMlF TypMFF Age2   Age3   Age4  
## TypeClear Valley -0.463                                          
## TypeMel Fry      -0.396  0.334                                   
## TypeMel Fry Free -0.463  0.721  0.334                            
## Age2             -0.533 -0.027  0.000 -0.027                     
## Age3             -0.319 -0.059  0.000 -0.059  0.272              
## Age4             -0.420 -0.081  0.000 -0.081  0.360  0.232       
## Age5             -0.380 -0.117  0.000 -0.117  0.341  0.223  0.295
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.57521136 -0.46418553  0.04419502  0.57764972  2.63882103 
## 
## Number of Observations: 1890
## Number of Groups: 945

# The above result shows the factor "Type" does not have a significant impact on the taste of French fries, but the factor "Age" does. This time every age level was significant. It appears taste improves significantly with age.
# Use diagnostic plots to check heteroscedasticity, normality, and influential observerations. These plots appears to indicate the model is a good fit.
plot(fit.Taste2); qqnorm(resid(fit.Taste2)); qqline(resid(fit.Taste2)) # diagnostic plots

anova(fit.Taste, fit.Taste2) # again (p=0.52) leaving out interactions was ok

##            Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## fit.Taste      1 22 6492.408 6614.383 -3224.204                        
## fit.Taste2     2 10 6479.513 6534.956 -3229.756 1 vs 2 11.10449    0.52

# make some plots
boxplot(Taste ~ Type, data) # again we don't see much difference among different types of oil

boxplot(Taste ~ Age, data) # apparently we can see Age 5 is significantly higher than other ages

# make a plot showing how taste change with ages of oil grouped by different types. We can see generally there is an uptrend with the increase of oil ages within each type of oil.
f <- ggplot(data=data.1, aes(x=Group.2, y=Taste)) 
f <- f + geom_bar(fill="grey", color="grey",stat = "identity") +
  guides(fill=F) + 
  ggtitle('Mean of Taste in Different Types and Ages of Oils') + 
  ylab('Taste') + 
  xlab('Age') + 
  theme(legend.position='none') + 
  facet_grid(. ~ Group.1) + 
  geom_text(aes(label=round(Liking, 0), size=1, hjust=0.5, vjust=-1))
print(f)

########### Texture ############
fit.Texture <- lme(Texture ~ Type*Age, data=data, random = ~ 1 | Subject, method='ML'); summary(fit.Texture)

## Linear mixed-effects model fit by maximum likelihood
##  Data: data 
##       AIC      BIC    logLik
##   6830.89 6952.865 -3393.445
## 
## Random effects:
##  Formula: ~1 | Subject
##         (Intercept) Residual
## StdDev:    1.010557 1.155458
## 
## Fixed effects: Texture ~ Type * Age 
##                           Value Std.Error  DF  t-value p-value
## (Intercept)            6.246512 0.1052464 940 59.35134  0.0000
## TypeClear Valley      -0.293387 0.1722857 930 -1.70291  0.0889
## TypeMel Fry           -0.051163 0.1120366 930 -0.45666  0.6480
## TypeMel Fry Free      -0.082449 0.1722857 930 -0.47856  0.6324
## Age2                  -0.086512 0.1641753 940 -0.52695  0.5984
## Age3                  -0.160797 0.2812827 940 -0.57166  0.5677
## Age4                  -0.010401 0.2101269 940 -0.04950  0.9605
## Age5                   0.696885 0.2366664 940  2.94459  0.0033
## TypeClear Valley:Age2  0.257196 0.2612281 930  0.98457  0.3251
## TypeMel Fry:Age2       0.137829 0.1747676 930  0.78864  0.4305
## TypeMel Fry Free:Age2  0.284354 0.2612281 930  1.08853  0.2766
## TypeClear Valley:Age3 -0.097883 0.4048192 930 -0.24179  0.8090
## TypeMel Fry:Age3       0.194020 0.2994305 930  0.64796  0.5172
## TypeMel Fry Free:Age3  0.024513 0.4048192 930  0.06055  0.9517
## TypeClear Valley:Age4  0.596749 0.3067484 930  1.94540  0.0520
## TypeMel Fry:Age4       0.162274 0.2236839 930  0.72546  0.4684
## TypeMel Fry Free:Age4  0.122654 0.3067484 930  0.39985  0.6894
## TypeClear Valley:Age5  0.056657 0.3261441 930  0.17372  0.8621
## TypeMel Fry:Age5       0.070031 0.2519356 930  0.27797  0.7811
## TypeMel Fry Free:Age5 -0.234280 0.3261441 930 -0.71833  0.4727
##  Correlation: 
##                       (Intr) TypClV TypMlF TypMFF Age2   Age3   Age4  
## TypeClear Valley      -0.611                                          
## TypeMel Fry           -0.532  0.325                                   
## TypeMel Fry Free      -0.611  0.645  0.325                            
## Age2                  -0.641  0.392  0.341  0.392                     
## Age3                  -0.374  0.229  0.199  0.229  0.240              
## Age4                  -0.501  0.306  0.267  0.306  0.321  0.187       
## Age5                  -0.445  0.272  0.237  0.272  0.285  0.166  0.223
## TypeClear Valley:Age2  0.403 -0.660 -0.214 -0.425 -0.628 -0.151 -0.202
## TypeMel Fry:Age2       0.341 -0.208 -0.641 -0.208 -0.532 -0.128 -0.171
## TypeMel Fry Free:Age2  0.403 -0.425 -0.214 -0.660 -0.628 -0.151 -0.202
## TypeClear Valley:Age3  0.260 -0.426 -0.138 -0.274 -0.167 -0.695 -0.130
## TypeMel Fry:Age3       0.199 -0.122 -0.374 -0.122 -0.128 -0.532 -0.100
## TypeMel Fry Free:Age3  0.260 -0.274 -0.138 -0.426 -0.167 -0.695 -0.130
## TypeClear Valley:Age4  0.343 -0.562 -0.183 -0.362 -0.220 -0.128 -0.685
## TypeMel Fry:Age4       0.267 -0.163 -0.501 -0.163 -0.171 -0.100 -0.532
## TypeMel Fry Free:Age4  0.343 -0.362 -0.183 -0.562 -0.220 -0.128 -0.685
## TypeClear Valley:Age5  0.323 -0.528 -0.172 -0.341 -0.207 -0.121 -0.162
## TypeMel Fry:Age5       0.237 -0.145 -0.445 -0.145 -0.152 -0.089 -0.119
## TypeMel Fry Free:Age5  0.323 -0.341 -0.172 -0.528 -0.207 -0.121 -0.162
##                       Age5   TCV:A2 TMF:A2 TMFF:A2 TCV:A3 TMF:A3 TMFF:A3
## TypeClear Valley                                                        
## TypeMel Fry                                                             
## TypeMel Fry Free                                                        
## Age2                                                                    
## Age3                                                                    
## Age4                                                                    
## Age5                                                                    
## TypeClear Valley:Age2 -0.179                                            
## TypeMel Fry:Age2      -0.152  0.335                                     
## TypeMel Fry Free:Age2 -0.179  0.657  0.335                              
## TypeClear Valley:Age3 -0.116  0.281  0.089  0.181                       
## TypeMel Fry:Age3      -0.089  0.080  0.240  0.080   0.370               
## TypeMel Fry Free:Age3 -0.116  0.181  0.089  0.281   0.707  0.370        
## TypeClear Valley:Age4 -0.153  0.370  0.117  0.239   0.239  0.068  0.154 
## TypeMel Fry:Age4      -0.119  0.107  0.321  0.107   0.069  0.187  0.069 
## TypeMel Fry Free:Age4 -0.153  0.239  0.117  0.370   0.154  0.068  0.239 
## TypeClear Valley:Age5 -0.726  0.348  0.110  0.225   0.225  0.064  0.145 
## TypeMel Fry:Age5      -0.532  0.095  0.285  0.095   0.062  0.166  0.062 
## TypeMel Fry Free:Age5 -0.726  0.225  0.110  0.348   0.145  0.064  0.225 
##                       TCV:A4 TMF:A4 TMFF:A4 TCV:A5 TMF:A5
## TypeClear Valley                                         
## TypeMel Fry                                              
## TypeMel Fry Free                                         
## Age2                                                     
## Age3                                                     
## Age4                                                     
## Age5                                                     
## TypeClear Valley:Age2                                    
## TypeMel Fry:Age2                                         
## TypeMel Fry Free:Age2                                    
## TypeClear Valley:Age3                                    
## TypeMel Fry:Age3                                         
## TypeMel Fry Free:Age3                                    
## TypeClear Valley:Age4                                    
## TypeMel Fry:Age4       0.365                             
## TypeMel Fry Free:Age4  0.699  0.365                      
## TypeClear Valley:Age5  0.297  0.086  0.191               
## TypeMel Fry:Age5       0.081  0.223  0.081   0.386       
## TypeMel Fry Free:Age5  0.191  0.086  0.297   0.732  0.386
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.59569377 -0.47422190  0.08868374  0.56862469  2.41823021 
## 
## Number of Observations: 1890
## Number of Groups: 945

# The p-value for the interaction term of 0.893 suggests we do not have to worry about interaction. So we can repeat the process but with a simple additive model.
fit.Texture2 <- lme(Texture ~ Type + Age, data=data, random = ~ 1 | Subject, method='ML'); summary(fit.Texture2)

## Linear mixed-effects model fit by maximum likelihood
##  Data: data 
##        AIC      BIC    logLik
##   6817.582 6873.025 -3398.791
## 
## Random effects:
##  Formula: ~1 | Subject
##         (Intercept) Residual
## StdDev:    1.010444 1.160205
## 
## Fixed effects: Texture ~ Type + Age 
##                      Value  Std.Error  DF  t-value p-value
## (Intercept)       6.178578 0.08536655 942 72.37704  0.0000
## TypeClear Valley -0.119680 0.10181227 942 -1.17549  0.2401
## TypeMel Fry       0.030476 0.07176145 942  0.42469  0.6712
## TypeMel Fry Free -0.029204 0.10181227 942 -0.28684  0.7743
## Age2              0.062718 0.10790045 940  0.58126  0.5612
## Age3             -0.141303 0.17045464 940 -0.82898  0.4073
## Age4              0.203587 0.12884799 940  1.58005  0.1144
## Age5              0.643887 0.13634103 940  4.72262  0.0000
##  Correlation: 
##                  (Intr) TypClV TypMlF TypMFF Age2   Age3   Age4  
## TypeClear Valley -0.469                                          
## TypeMel Fry      -0.420  0.352                                   
## TypeMel Fry Free -0.469  0.689  0.352                            
## Age2             -0.527 -0.026  0.000 -0.026                     
## Age3             -0.315 -0.058  0.000 -0.058  0.272              
## Age4             -0.415 -0.080  0.000 -0.080  0.360  0.232       
## Age5             -0.375 -0.114  0.000 -0.114  0.341  0.223  0.295
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.60986958 -0.48680611  0.09420918  0.56192462  2.44981419 
## 
## Number of Observations: 1890
## Number of Groups: 945

# The above result shows the factor "Type" does not have a significant impact on the Texture of French fries, but the factor "Age" does. Only age 5 was significant compared with age 1.
# Use diagnostic plots to check heteroscedasticity, normality, and influential observerations. These plots appears to indicate the model is a good fit.
plot(fit.Texture2); qqnorm(resid(fit.Texture2)); qqline(resid(fit.Texture2)) # diagnostic plots

anova(fit.Texture, fit.Texture2) # again (p=0.55) leaving interaction terms out was ok

##              Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## fit.Texture      1 22 6830.890 6952.865 -3393.445                        
## fit.Texture2     2 10 6817.582 6873.025 -3398.791 1 vs 2 10.69158  0.5555

# make some plots
boxplot(Texture ~ Type, data) # again we don't see much difference among different types of oil

boxplot(Texture ~ Age, data) # apparently we can see Age 5 is significantly higher than other ages

# make a plot showing how Texture change with ages of oil grouped by different types. We can see generally there is an uptrend with the increase of oil ages within each type of oil.
g <- ggplot(data=data.1, aes(x=Group.2, y=Texture)) 
g <- g + geom_bar(fill="grey", color="grey",stat = "identity") +
  guides(fill=F) + 
  ggtitle('Mean of Texture in Different Types and Ages of Oils') + 
  ylab('Texture') + 
  xlab('Age') + 
  theme(legend.position='none') + 
  facet_grid(. ~ Group.1) + 
  geom_text(aes(label=round(Liking, 0), size=1, hjust=0.5, vjust=-1))
print(g)

########## Color ##########
# residuals were pretty non-Gaussian, so I used a log transform on the response. This helped, but they're still not great. Fortunately we have a large sample size, so I think we might still be able to use the model.
fit.Color <- lme(log(Color) ~ Type*Age, data=data, random = ~ 1 | Subject, method='ML'); summary(fit.Color)

## Linear mixed-effects model fit by maximum likelihood
##  Data: data 
##         AIC       BIC   logLik
##   -799.2029 -677.2276 421.6015
## 
## Random effects:
##  Formula: ~1 | Subject
##         (Intercept)  Residual
## StdDev:  0.09860937 0.1702778
## 
## Fixed effects: log(Color) ~ Type * Age 
##                            Value  Std.Error  DF   t-value p-value
## (Intercept)            1.4956272 0.01349116 940 110.85979  0.0000
## TypeClear Valley      -0.0267115 0.02208469 930  -1.20950  0.2268
## TypeMel Fry           -0.0068329 0.01651064 930  -0.41385  0.6791
## TypeMel Fry Free       0.0021961 0.02208469 930   0.09944  0.9208
## Age2                   0.0167624 0.02104506 940   0.79650  0.4259
## Age3                  -0.0225970 0.03605664 940  -0.62671  0.5310
## Age4                   0.0245811 0.02693544 940   0.91259  0.3617
## Age5                   0.0144147 0.03033743 940   0.47514  0.6348
## TypeClear Valley:Age2  0.0203515 0.03348592 930   0.60776  0.5435
## TypeMel Fry:Age2       0.0114539 0.02575519 930   0.44472  0.6566
## TypeMel Fry Free:Age2  0.0022957 0.03348592 930   0.06856  0.9454
## TypeClear Valley:Age3  0.0345926 0.05189235 930   0.66662  0.5052
## TypeMel Fry:Age3      -0.0339419 0.04412655 930  -0.76919  0.4420
## TypeMel Fry Free:Age3  0.0154689 0.05189235 930   0.29810  0.7657
## TypeClear Valley:Age4 -0.0087719 0.03932100 930  -0.22309  0.8235
## TypeMel Fry:Age4       0.0291862 0.03296391 930   0.88540  0.3762
## TypeMel Fry Free:Age4 -0.0144923 0.03932100 930  -0.36856  0.7125
## TypeClear Valley:Age5 -0.0298809 0.04180727 930  -0.71473  0.4750
## TypeMel Fry:Age5       0.0426276 0.03712731 930   1.14815  0.2512
## TypeMel Fry Free:Age5 -0.0584722 0.04180727 930  -1.39861  0.1623
##  Correlation: 
##                       (Intr) TypClV TypMlF TypMFF Age2   Age3   Age4  
## TypeClear Valley      -0.611                                          
## TypeMel Fry           -0.612  0.374                                   
## TypeMel Fry Free      -0.611  0.531  0.374                            
## Age2                  -0.641  0.392  0.392  0.392                     
## Age3                  -0.374  0.229  0.229  0.229  0.240              
## Age4                  -0.501  0.306  0.306  0.306  0.321  0.187       
## Age5                  -0.445  0.272  0.272  0.272  0.285  0.166  0.223
## TypeClear Valley:Age2  0.403 -0.660 -0.247 -0.350 -0.628 -0.151 -0.202
## TypeMel Fry:Age2       0.392 -0.240 -0.641 -0.240 -0.612 -0.147 -0.196
## TypeMel Fry Free:Age2  0.403 -0.350 -0.247 -0.660 -0.628 -0.151 -0.202
## TypeClear Valley:Age3  0.260 -0.426 -0.159 -0.226 -0.167 -0.695 -0.130
## TypeMel Fry:Age3       0.229 -0.140 -0.374 -0.140 -0.147 -0.612 -0.115
## TypeMel Fry Free:Age3  0.260 -0.226 -0.159 -0.426 -0.167 -0.695 -0.130
## TypeClear Valley:Age4  0.343 -0.562 -0.210 -0.298 -0.220 -0.128 -0.685
## TypeMel Fry:Age4       0.306 -0.187 -0.501 -0.187 -0.196 -0.115 -0.612
## TypeMel Fry Free:Age4  0.343 -0.298 -0.210 -0.562 -0.220 -0.128 -0.685
## TypeClear Valley:Age5  0.323 -0.528 -0.197 -0.280 -0.207 -0.121 -0.162
## TypeMel Fry:Age5       0.272 -0.166 -0.445 -0.166 -0.174 -0.102 -0.136
## TypeMel Fry Free:Age5  0.323 -0.280 -0.197 -0.528 -0.207 -0.121 -0.162
##                       Age5   TCV:A2 TMF:A2 TMFF:A2 TCV:A3 TMF:A3 TMFF:A3
## TypeClear Valley                                                        
## TypeMel Fry                                                             
## TypeMel Fry Free                                                        
## Age2                                                                    
## Age3                                                                    
## Age4                                                                    
## Age5                                                                    
## TypeClear Valley:Age2 -0.179                                            
## TypeMel Fry:Age2      -0.174  0.385                                     
## TypeMel Fry Free:Age2 -0.179  0.547  0.385                              
## TypeClear Valley:Age3 -0.116  0.281  0.102  0.149                       
## TypeMel Fry:Age3      -0.102  0.092  0.240  0.092   0.425               
## TypeMel Fry Free:Age3 -0.116  0.149  0.102  0.281   0.613  0.425        
## TypeClear Valley:Age4 -0.153  0.370  0.135  0.197   0.239  0.079  0.127 
## TypeMel Fry:Age4      -0.136  0.123  0.321  0.123   0.080  0.187  0.080 
## TypeMel Fry Free:Age4 -0.153  0.197  0.135  0.370   0.127  0.079  0.239 
## TypeClear Valley:Age5 -0.726  0.348  0.127  0.185   0.225  0.074  0.119 
## TypeMel Fry:Age5      -0.612  0.110  0.285  0.110   0.071  0.166  0.071 
## TypeMel Fry Free:Age5 -0.726  0.185  0.127  0.348   0.119  0.074  0.225 
##                       TCV:A4 TMF:A4 TMFF:A4 TCV:A5 TMF:A5
## TypeClear Valley                                         
## TypeMel Fry                                              
## TypeMel Fry Free                                         
## Age2                                                     
## Age3                                                     
## Age4                                                     
## Age5                                                     
## TypeClear Valley:Age2                                    
## TypeMel Fry:Age2                                         
## TypeMel Fry Free:Age2                                    
## TypeClear Valley:Age3                                    
## TypeMel Fry:Age3                                         
## TypeMel Fry Free:Age3                                    
## TypeClear Valley:Age4                                    
## TypeMel Fry:Age4       0.419                             
## TypeMel Fry Free:Age4  0.603  0.419                      
## TypeClear Valley:Age5  0.297  0.099  0.157               
## TypeMel Fry:Age5       0.093  0.223  0.093   0.444       
## TypeMel Fry Free:Age5  0.157  0.099  0.297   0.645  0.444
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -7.0265412 -0.4261101  0.3399997  0.4657097  2.4459530 
## 
## Number of Observations: 1890
## Number of Groups: 945

# No significant interaction terms suggests we do not have to worry about interaction. So we can repeat the process but with a simple additive model.
fit.Color2 <- lme(log(Color) ~ Type + Age, data=data, random = ~ 1 | Subject, method='ML'); summary(fit.Color2)

## Linear mixed-effects model fit by maximum likelihood
##  Data: data 
##         AIC       BIC   logLik
##   -812.1895 -756.7462 416.0948
## 
## Random effects:
##  Formula: ~1 | Subject
##         (Intercept) Residual
## StdDev:  0.09925574 0.170619
## 
## Fixed effects: log(Color) ~ Type + Age 
##                       Value  Std.Error  DF   t-value p-value
## (Intercept)       1.4941025 0.01069060 942 139.75856  0.0000
## TypeClear Valley -0.0243026 0.01304825 942  -1.86252  0.0628
## TypeMel Fry       0.0024829 0.01055319 942   0.23527  0.8140
## TypeMel Fry Free -0.0076888 0.01304825 942  -0.58926  0.5558
## Age2              0.0251178 0.01295151 940   1.93937  0.0528
## Age3             -0.0171474 0.02046001 940  -0.83809  0.4022
## Age4              0.0268855 0.01546588 940   1.73837  0.0825
## Age5             -0.0008585 0.01636529 940  -0.05246  0.9582
##  Correlation: 
##                  (Intr) TypClV TypMlF TypMFF Age2   Age3   Age4  
## TypeClear Valley -0.487                                          
## TypeMel Fry      -0.494  0.404                                   
## TypeMel Fry Free -0.487  0.591  0.404                            
## Age2             -0.505 -0.025  0.000 -0.025                     
## Age3             -0.302 -0.054  0.000 -0.054  0.272              
## Age4             -0.398 -0.075  0.000 -0.075  0.360  0.232       
## Age5             -0.360 -0.107  0.000 -0.107  0.341  0.223  0.295
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -7.1158267 -0.3884339  0.3182720  0.4199311  2.4144825 
## 
## Number of Observations: 1890
## Number of Groups: 945

# The above result shows neither type nor age has any significant levels.
# Use diagnostic plots to check heteroscedasticity, normality, and influential observerations. These plots aren't great.
plot(fit.Color2); qqnorm(resid(fit.Color2)); qqline(resid(fit.Color2)) # diagnostic plots (note the qqplot isn't very good)

anova(fit.Color, fit.Color2) # again (p=0.52) leaving out interactions was ok

##            Model df       AIC       BIC   logLik   Test  L.Ratio p-value
## fit.Color      1 22 -799.2029 -677.2276 421.6015                        
## fit.Color2     2 10 -812.1895 -756.7462 416.0948 1 vs 2 11.01341  0.5278

# make some plots
boxplot(Color ~ Type, data) # again we don't see much difference among different types of oil

boxplot(Color ~ Age, data) # not much difference here either

# make a plot showing how color changes with ages of oil grouped by different types. There isn't much of a trend at all. I'm not sure we can conclude anything about color.
h <- ggplot(data=data.1, aes(x=Group.2, y=Color)) 
h <- h + geom_bar(fill="grey", color="grey",stat = "identity") +
  guides(fill=F) + 
  ggtitle('Mean of Color in Different Types and Ages of Oils') + 
  ylab('Color') + 
  xlab('Age') + 
  theme(legend.position='none') + 
  facet_grid(. ~ Group.1) + 
  geom_text(aes(label=round(Liking, 0), size=1, hjust=0.5, vjust=-1))
print(h)

########## Temperature ##########
# not sure we can trust this model either
fit.Temperature <- lme(Temperature ~ Type*Age, data=data, random = ~ 1 | Subject, method='ML'); summary(fit.Temperature)

## Linear mixed-effects model fit by maximum likelihood
##  Data: data 
##       AIC      BIC   logLik
##   3804.58 3926.556 -1880.29
## 
## Random effects:
##  Formula: ~1 | Subject
##         (Intercept)  Residual
## StdDev:   0.4050867 0.5426561
## 
## Fixed effects: Temperature ~ Type * Age 
##                           Value  Std.Error  DF  t-value p-value
## (Intercept)            4.446512 0.04642951 940 95.76909  0.0000
## TypeClear Valley      -0.149637 0.07600395 930 -1.96880  0.0493
## TypeMel Fry            0.037209 0.05261754 930  0.70717  0.4796
## TypeMel Fry Free      -0.040262 0.07600395 930 -0.52973  0.5964
## Age2                  -0.073178 0.07242607 940 -1.01039  0.3126
## Age3                  -0.103654 0.12408808 940 -0.83533  0.4037
## Age4                   0.164599 0.09269766 940  1.77566  0.0761
## Age5                   0.176130 0.10440554 940  1.68698  0.0919
## TypeClear Valley:Age2  0.338208 0.11524100 930  2.93479  0.0034
## TypeMel Fry:Age2       0.142791 0.08207887 930  1.73968  0.0822
## TypeMel Fry Free:Age2  0.285976 0.11524100 930  2.48155  0.0133
## TypeClear Valley:Age3  0.362335 0.17858631 930  2.02891  0.0428
## TypeMel Fry:Age3       0.048505 0.14062629 930  0.34492  0.7302
## TypeMel Fry Free:Age3  0.252960 0.17858631 930  1.41646  0.1570
## TypeClear Valley:Age4  0.275368 0.13532227 930  2.03490  0.0421
## TypeMel Fry:Age4       0.060013 0.10505222 930  0.57127  0.5680
## TypeMel Fry Free:Age4  0.231782 0.13532227 930  1.71282  0.0871
## TypeClear Valley:Age5  0.273662 0.14387873 930  1.90203  0.0575
## TypeMel Fry:Age5       0.132602 0.11832050 930  1.12070  0.2627
## TypeMel Fry Free:Age5  0.310953 0.14387873 930  2.16122  0.0309
##  Correlation: 
##                       (Intr) TypClV TypMlF TypMFF Age2   Age3   Age4  
## TypeClear Valley      -0.611                                          
## TypeMel Fry           -0.567  0.346                                   
## TypeMel Fry Free      -0.611  0.597  0.346                            
## Age2                  -0.641  0.392  0.363  0.392                     
## Age3                  -0.374  0.229  0.212  0.229  0.240              
## Age4                  -0.501  0.306  0.284  0.306  0.321  0.187       
## Age5                  -0.445  0.272  0.252  0.272  0.285  0.166  0.223
## TypeClear Valley:Age2  0.403 -0.660 -0.228 -0.394 -0.628 -0.151 -0.202
## TypeMel Fry:Age2       0.363 -0.222 -0.641 -0.222 -0.567 -0.136 -0.182
## TypeMel Fry Free:Age2  0.403 -0.394 -0.228 -0.660 -0.628 -0.151 -0.202
## TypeClear Valley:Age3  0.260 -0.426 -0.147 -0.254 -0.167 -0.695 -0.130
## TypeMel Fry:Age3       0.212 -0.130 -0.374 -0.130 -0.136 -0.567 -0.106
## TypeMel Fry Free:Age3  0.260 -0.254 -0.147 -0.426 -0.167 -0.695 -0.130
## TypeClear Valley:Age4  0.343 -0.562 -0.194 -0.336 -0.220 -0.128 -0.685
## TypeMel Fry:Age4       0.284 -0.173 -0.501 -0.173 -0.182 -0.106 -0.567
## TypeMel Fry Free:Age4  0.343 -0.336 -0.194 -0.562 -0.220 -0.128 -0.685
## TypeClear Valley:Age5  0.323 -0.528 -0.183 -0.316 -0.207 -0.121 -0.162
## TypeMel Fry:Age5       0.252 -0.154 -0.445 -0.154 -0.162 -0.094 -0.126
## TypeMel Fry Free:Age5  0.323 -0.316 -0.183 -0.528 -0.207 -0.121 -0.162
##                       Age5   TCV:A2 TMF:A2 TMFF:A2 TCV:A3 TMF:A3 TMFF:A3
## TypeClear Valley                                                        
## TypeMel Fry                                                             
## TypeMel Fry Free                                                        
## Age2                                                                    
## Age3                                                                    
## Age4                                                                    
## Age5                                                                    
## TypeClear Valley:Age2 -0.179                                            
## TypeMel Fry:Age2      -0.162  0.356                                     
## TypeMel Fry Free:Age2 -0.179  0.611  0.356                              
## TypeClear Valley:Age3 -0.116  0.281  0.094  0.168                       
## TypeMel Fry:Age3      -0.094  0.085  0.240  0.085   0.394               
## TypeMel Fry Free:Age3 -0.116  0.168  0.094  0.281   0.668  0.394        
## TypeClear Valley:Age4 -0.153  0.370  0.125  0.221   0.239  0.073  0.143 
## TypeMel Fry:Age4      -0.126  0.114  0.321  0.114   0.074  0.187  0.074 
## TypeMel Fry Free:Age4 -0.153  0.221  0.125  0.370   0.143  0.073  0.239 
## TypeClear Valley:Age5 -0.726  0.348  0.117  0.208   0.225  0.068  0.134 
## TypeMel Fry:Age5      -0.567  0.102  0.285  0.102   0.066  0.166  0.066 
## TypeMel Fry Free:Age5 -0.726  0.208  0.117  0.348   0.134  0.068  0.225 
##                       TCV:A4 TMF:A4 TMFF:A4 TCV:A5 TMF:A5
## TypeClear Valley                                         
## TypeMel Fry                                              
## TypeMel Fry Free                                         
## Age2                                                     
## Age3                                                     
## Age4                                                     
## Age5                                                     
## TypeClear Valley:Age2                                    
## TypeMel Fry:Age2                                         
## TypeMel Fry Free:Age2                                    
## TypeClear Valley:Age3                                    
## TypeMel Fry:Age3                                         
## TypeMel Fry Free:Age3                                    
## TypeClear Valley:Age4                                    
## TypeMel Fry:Age4       0.388                             
## TypeMel Fry Free:Age4  0.659  0.388                      
## TypeClear Valley:Age5  0.297  0.092  0.177               
## TypeMel Fry:Age5       0.086  0.223  0.086   0.411       
## TypeMel Fry Free:Age5  0.177  0.092  0.297   0.696  0.411
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -4.5556069 -0.4841690  0.2920085  0.5004378  2.0712546 
## 
## Number of Observations: 1890
## Number of Groups: 945

# No significant interaction terms suggests we do not have to worry about interaction. So we can repeat the process but with a simple additive model.
fit.Temperature2 <- lme(Temperature ~ Type + Age, data=data, random = ~ 1 | Subject, method='ML'); summary(fit.Temperature2)

## Linear mixed-effects model fit by maximum likelihood
##  Data: data 
##        AIC      BIC    logLik
##   3794.923 3850.366 -1887.462
## 
## Random effects:
##  Formula: ~1 | Subject
##         (Intercept)  Residual
## StdDev:    0.408085 0.5440194
## 
## Fixed effects: Temperature ~ Type + Age 
##                     Value  Std.Error  DF   t-value p-value
## (Intercept)      4.356049 0.03737669 942 116.54453  0.0000
## TypeClear Valley 0.049472 0.04498420 942   1.09976  0.2717
## TypeMel Fry      0.102857 0.03364890 942   3.05678  0.0023
## TypeMel Fry Free 0.135186 0.04498420 942   3.00519  0.0027
## Age2             0.091368 0.04649226 940   1.96522  0.0497
## Age3             0.043615 0.07344567 940   0.59384  0.5528
## Age4             0.287736 0.05551816 940   5.18273  0.0000
## Age5             0.341995 0.05874677 940   5.82151  0.0000
##  Correlation: 
##                  (Intr) TypClV TypMlF TypMFF Age2   Age3   Age4  
## TypeClear Valley -0.476                                          
## TypeMel Fry      -0.450  0.374                                   
## TypeMel Fry Free -0.476  0.650  0.374                            
## Age2             -0.518 -0.026  0.000 -0.026                     
## Age3             -0.310 -0.056  0.000 -0.056  0.272              
## Age4             -0.409 -0.078  0.000 -0.078  0.360  0.232       
## Age5             -0.369 -0.112  0.000 -0.112  0.341  0.223  0.295
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -4.5588341 -0.5061186  0.3070041  0.5279649  2.0669454 
## 
## Number of Observations: 1890
## Number of Groups: 945

# The above result shows the Mel Fry types had an impact on the temperature of French fries, as well as ages 2, 4, and 5. Strange result
# Use diagnostic plots to check heteroscedasticity, normality, and influential observerations.
plot(fit.Temperature2); qqnorm(resid(fit.Temperature2)); qqline(resid(fit.Temperature2)) # diagnostic plots. Residuals don't look good.

anova(fit.Temperature, fit.Temperature2) # again (p=0.27) leaving out interactions was ok

##                  Model df      AIC      BIC    logLik   Test  L.Ratio
## fit.Temperature      1 22 3804.580 3926.556 -1880.290                
## fit.Temperature2     2 10 3794.923 3850.366 -1887.461 1 vs 2 14.34262
##                  p-value
## fit.Temperature         
## fit.Temperature2  0.2794

# make some plots
boxplot(Temperature ~ Type, data) # again we don't see much difference among different types of oil

boxplot(Temperature ~ Age, data) # does look like liking increases with temperature, contrary to model

# make a plot showing how Temperaure change with ages of oil grouped by different types. We can see generally there is an uptrend with the increase of oil ages within each type of oil.
j <- ggplot(data=data.1, aes(x=Group.2, y=Temperature)) 
j <- j + geom_bar(fill="grey", color="grey",stat = "identity") +
  guides(fill=F) + 
  ggtitle('Mean of Temperature in Different Types and Ages of Oils') + 
  ylab('Temperature') + 
  xlab('Age') + 
  theme(legend.position='none') + 
  facet_grid(. ~ Group.1) + 
  geom_text(aes(label=round(Liking, 0), size=1, hjust=0.5, vjust=-1))
print(j)

Question 3: How are the characteristics (temperature, appearance, color, taste and texture) of French fries related to the overall liking?

# Model 1: Least squares regression model 
data.2 <- data[, 3:9]
OLS <- lme(Liking~Appearance+Color+Taste+Texture+Temperature, data=data.2, random = ~ 1 | Subject, method='ML'); summary(OLS)

## Linear mixed-effects model fit by maximum likelihood
##  Data: data.2 
##       AIC      BIC    logLik
##   3344.31 3388.665 -1664.155
## 
## Random effects:
##  Formula: ~1 | Subject
##         (Intercept)  Residual
## StdDev:   0.2798588 0.5207953
## 
## Fixed effects: Liking ~ Appearance + Color + Taste + Texture + Temperature 
##                  Value  Std.Error  DF  t-value p-value
## (Intercept) -0.3835078 0.13245597 944 -2.89536  0.0039
## Appearance   0.1353313 0.01307124 940 10.35336  0.0000
## Color        0.0615116 0.02020927 940  3.04373  0.0024
## Taste        0.5575995 0.01367111 940 40.78669  0.0000
## Texture      0.2687682 0.01291434 940 20.81161  0.0000
## Temperature  0.0841870 0.02023759 940  4.15993  0.0000
##  Correlation: 
##             (Intr) Apprnc Color  Taste  Textur
## Appearance  -0.179                            
## Color       -0.540 -0.193                     
## Taste       -0.089 -0.278  0.016              
## Texture     -0.030 -0.297 -0.023 -0.530       
## Temperature -0.610  0.051 -0.038 -0.087 -0.056
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -4.36706334 -0.43255104 -0.03315379  0.43280105  4.38939190 
## 
## Number of Observations: 1890
## Number of Groups: 945

plot(OLS); qqnorm(resid(OLS)); qqline(resid(OLS)) # not good (I actually think these look pretty okay - John)

library(car)
vif(OLS) # all ok - multicollinearity not an issue

##  Appearance       Color       Taste     Texture Temperature 
##    1.646077    1.067470    2.031741    2.060996    1.027350

# Model 2: Best subset selection using BIC to select the number of predictors
library(leaps)
BIC <- regsubsets(Liking~., data.2, nvmax=5) 
summary(BIC)

## Subset selection object
## Call: regsubsets.formula(Liking ~ ., data.2, nvmax = 5)
## 6 Variables  (and intercept)
##             Forced in Forced out
## Temperature     FALSE      FALSE
## Appearance      FALSE      FALSE
## Color           FALSE      FALSE
## Taste           FALSE      FALSE
## Texture         FALSE      FALSE
## Subject         FALSE      FALSE
## 1 subsets of each size up to 5
## Selection Algorithm: exhaustive
##          Temperature Appearance Color Taste Texture Subject
## 1  ( 1 ) " "         " "        " "   "*"   " "     " "    
## 2  ( 1 ) " "         " "        " "   "*"   "*"     " "    
## 3  ( 1 ) " "         "*"        " "   "*"   "*"     " "    
## 4  ( 1 ) "*"         "*"        " "   "*"   "*"     " "    
## 5  ( 1 ) "*"         "*"        "*"   "*"   "*"     " "

which.min(summary(BIC)$bic) # determine the number of variables with the lowest BIC value; 5 is result

## [1] 5

coef(BIC, id=5) 

## (Intercept) Temperature  Appearance       Color       Taste     Texture 
## -0.36939032  0.07833715  0.13406545  0.06021952  0.56229244  0.26836336

# Model 3: Stepwise Regression
library(MASS)

## 
## Attaching package: 'MASS'
## 
## The following object is masked from 'package:dplyr':
## 
##     select

fit.step <- lm(Liking~., data.2)
step <- stepAIC(fit.step, direction="both")

## Start:  AIC=-1973
## Liking ~ Temperature + Appearance + Color + Taste + Texture + 
##     Subject
## 
##               Df Sum of Sq     RSS      AIC
## - Subject      1      0.06  660.57 -1974.83
## <none>                      660.51 -1973.00
## - Color        1      3.20  663.71 -1965.85
## - Temperature  1      5.37  665.88 -1959.69
## - Appearance   1     37.56  698.07 -1870.46
## - Texture      1    152.49  812.99 -1582.42
## - Taste        1    606.41 1266.92  -743.97
## 
## Step:  AIC=-1974.83
## Liking ~ Temperature + Appearance + Color + Taste + Texture
## 
##               Df Sum of Sq     RSS      AIC
## <none>                      660.57 -1974.83
## + Subject      1      0.06  660.51 -1973.00
## - Color        1      3.17  663.73 -1967.80
## - Temperature  1      5.45  666.02 -1961.30
## - Appearance   1     37.70  698.27 -1871.93
## - Texture      1    152.50  813.06 -1584.26
## - Taste        1    606.69 1267.26  -745.47

step$anova # display results

## Stepwise Model Path 
## Analysis of Deviance Table
## 
## Initial Model:
## Liking ~ Temperature + Appearance + Color + Taste + Texture + 
##     Subject
## 
## Final Model:
## Liking ~ Temperature + Appearance + Color + Taste + Texture
## 
## 
##        Step Df   Deviance Resid. Df Resid. Dev       AIC
## 1                              1883   660.5091 -1972.997
## 2 - Subject  1 0.05715246      1884   660.5662 -1974.834

### Looks like each of the three models includes all five variables.